my_stocks <- c("ITUB3.SA", "ENBR3.SA",
"FLRY3.SA", "B3SA3.SA", "BBDC3.SA",
"EGIE3.SA", "TAEE11.SA", "WEGE3.SA",
"PSSA3.SA", "^BVSP")
all_asset_returns <-
BatchGetSymbols(my_stocks,
first.date = "2000-01-01",
bench.ticker = "^BVSP",
type.return = "log",
freq.data = "daily"
)$df.tickers %>%
select(ref.date, ticker, ret.closing.prices, price.close) %>%
drop_na()
tickers_filter_table = tibble(
stock = all_asset_returns$ticker %>%
unique,
add_date = c(0, 2300, 2000, 0, 3000, 4000, 0)
)
do_filter_ref_date <- function(stock, filter_table) {
stock %>%
filter(ref.date > (
stock$ref.date[1] +
filter_table$add_date[filter_table$stock == stock$ticker[1]]
)
)
}
asset_returns <- read_csv("\n", col_names = names(all_asset_returns)) %>%
mutate(
ref.date = as.Date(ref.date),
ret.closing.prices = as.numeric(),
price.close = as.numeric(),
)
for (t in unique(all_asset_returns$ticker)) {
asset_returns <- asset_returns %>%
add_row(all_asset_returns %>%
filter(ticker == t) %>%
drop_na('ret.closing.prices') %>%
do_filter_ref_date(tickers_filter_table))
}
plot_returns <- function(stock) {
if (nrow(stock) > 0) {
df_plot <- stock
par(mfrow=c(1, 3))
ts.plot(df_plot$price.close, main = paste("Ticker", stock$ticker[1], " closing prices"))
ts.plot(
df_plot$ret.closing.prices,
main = paste("Ticker", stock$ticker[1], "ret")
)
ts.plot(
df_plot$ret.closing.prices^2,
main = paste("Ticker", stock$ticker[1], "ret^2")
)
}
}
plot_asset_returns <- function (asset_returns) {
tickers <- asset_returns$ticker %>%
unique
for (t in tickers) {
asset_returns %>%
filter(ticker == t) %>%
plot_returns()
}
}
box_test <- function (df, lag = 1, type = 'Ljung-Box') Box.test(df, lag = lag, type = type)
descriptive_asset_evaluation <- function (asset_returns) {
tickers <- asset_returns$ticker %>%
unique
for (t in tickers) {
df_box <- asset_returns %>%
filter(ticker == t)
acf(df_box$ret.closing.prices, main = paste("ACF", t)) %>%
autoplot() %>%
theme(
plot.title = element_text(vjust = -50)
)
print(box_test(df_box$ret.closing.prices))
print(adf.test(df_box$ret.closing.prices))
print(ArchTest(df_box$ret.closing.prices))
}
}
visualize_residual_distribution <- function(asset_returns) {
tickers <- asset_returns$ticker %>%
unique
par(mfrow = c(1,2))
for (t in tickers) {
ret <- (asset_returns %>%
filter(ticker == t) %>%
drop_na('ret.closing.prices'))$ret.closing.prices
h <- hist(ret,
breaks=20,
col="red",
xlab="",
main= paste("Histogram", t)
)
xfit <- seq(min(ret),
max(ret),
length = 40)
yfit <- dnorm(xfit,
mean = mean(ret),
sd = sd(ret))
yfit <- yfit * diff(h$mids[1:2]) * length(ret)
lines(xfit, yfit, col="blue", lwd=2)
qqnorm(ret, pch = 1, frame = FALSE)
qqline(ret, col = "steelblue", lwd = 2)
print(paste("Realizando teste de normalidade para o ticker", t))
print(shapiro.test(as.vector(ret[c(1:4999)])))
}
}
calculate_assets_beta <- function (asset_returns) {
asset_returns %>%
group_by(ticker) %>%
do(model = lm(ret.closing.prices ~ market_ret$ret.closing.prices,
data = .))
}
plot_asset_returns(asset_returns)
Os retornos aparentemente são estacionários. Iremos fazer o teste de ljung-box para confirmar o ruído branco nos resíduos e também realizar o teste de raíz unitária.
descriptive_asset_evaluation(asset_returns)
##
## Box-Ljung test
##
## data: df
## X-squared = 8.4382, df = 1, p-value = 0.003674
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -17.489, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 273.49, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 40.974, df = 1, p-value = 1.543e-10
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -16.892, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 527.07, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 3.0717, df = 1, p-value = 0.07967
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -14.665, Lag order = 16, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 824.95, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 23.779, df = 1, p-value = 1.081e-06
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -17.779, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 387.87, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 3.2264, df = 1, p-value = 0.07246
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -14.784, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 642.15, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 2.8955, df = 1, p-value = 0.08883
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -10.607, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 128.7, df = 12, p-value < 2.2e-16
##
## Box-Ljung test
##
## data: df
## X-squared = 3.4689, df = 1, p-value = 0.06253
##
##
## Augmented Dickey-Fuller Test
##
## data: df_box$ret.closing.prices
## Dickey-Fuller = -16.147, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary
##
##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: df_box$ret.closing.prices
## Chi-squared = 1531.6, df = 12, p-value < 2.2e-16
visualize_residual_distribution(asset_returns)
calculate_portfolio_covariability <- function(tib_returns) {
cov(tib_returns$returns, tib_returns$ret.closing.prices) / var(tib_returns$ret.closing.prices)
}
fit_ts_stock <- function (stock) {
if (nrow(stock) > 0) {
stock <- stock %>%
timetk::future_frame(
.length_out = 3,
.date_var = ref.date,
.bind_data = TRUE
) %>%
mutate(ref.date = as.POSIXct(ref.date))
stock_train <- stock %>%
drop_na()
stock_future <- stock %>%
filter(is.na(ret.closing.prices))
stock_model <- garchmodels::garch_reg(
mode = "regression",
arch_order = tune::tune(),
garch_order = tune::tune(),
ma_order = tune::tune(),
ar_order = tune::tune(),
tune_by = "sigmaFor"
) %>%
parsnip::set_engine("rugarch")
stock_recipe <- recipes::recipe(
ret.closing.prices ~ ref.date,
data = stock_train
)
stock_wflw <- workflow() %>%
add_recipe(stock_recipe) %>%
add_model(stock_model)
stock_resamples <- time_series_cv(
stock_train,
date_var = ref.date,
initial = "1 year",
assess = "3 months",
skip = "1 month",
cumulative = TRUE,
slice_limit = 5
)
stock_tune_results <- tune_grid(
object = stock_wflw,
resamples = stock_resamples,
param_info = dials::parameters(stock_wflw),
grid = 5,
control = control_grid(
verbose = TRUE,
allow_par = TRUE,
parallel_over = "everything"
)
)
stock_tune_results %>%
tune::show_best(metric = "rmse")
}
}
results <- asset_returns %>%
group_split(ticker) %>%
map(fit_ts_stock)
results
evaluate_assets <- function(tickers,
bench_ticker = "^BVSP",
market = "^BVSP",
first.date = "2000-01-01",
last.date = Sys.Date(),
assets_weights = NULL
) {
asset_returns <- BatchGetSymbols(tickers,
first.date = first.date,
last.date = last.date,
bench.ticker = bench_ticker,
type.return = "log",
freq.data = "daily"
)$df.tickers %>%
select(ref.date, ticker, ret.closing.prices) %>%
drop_na()
tickers <- asset_returns$ticker %>%
unique
for (t in tickers) {
asset_returns %>%
filter(ticker == t) %>%
select(ref.date, ret.closing.prices) %>%
drop_na()
}
portfolio_ret <- asset_returns %>%
tq_portfolio(assets_col = ticker,
returns_col = ret.closing.prices,
weights = assets_weights,
col_rename = "returns")
market_ret <-
BatchGetSymbols(market,
first.date = first.date,
last.date = last.date,
freq.data = "daily",
type.return = "log")$df.tickers %>%
select(ref.date, ret.closing.prices) %>%
as_tibble()
ts.plot(portfolio_ret$returns)
ts.plot(market_ret$ret.closing.prices)
portfolio_market_ts <- portfolio_ret %>%
inner_join(market_ret, by = 'ref.date')
calculate_portfolio_covariability(portfolio_market_ts)
betas <- asset_returns %>%
group_by(ticker) %>%
inner_join(market_ret, by = 'ref.date') %>%
do(model = lm(ret.closing.prices.x ~ ret.closing.prices.y,
data = .))
betas
}
assets_lm <- evaluate_assets(my_stocks)
assets_lm$betas <- map(assets_lm$model, coef) %>%
map_dbl(2)
assets_lm
## para calcular o beta do portfolio precisamos dos weights de cada uma das ações, e então multiplicar cada beta com o seu peso respectivo - leia-se beta_porfolio_byhand do arquivo capm_model.rmd
dlm_routine <- function (df, ticker_string, index_string) {
my_dlm = function(parm, x.mat) {
parm = exp(parm)
return (
dlmModReg(
X = x.mat,
dV = parm[1],
dW = c(parm[2], parm[3])
)
)
}
ticker_ret <- df %>%
dplyr::select('ref.date','ticker','ret.closing.prices') %>%
dplyr::filter(ticker == ticker_string) %>%
drop_na('ret.closing.prices')
index_ret <- df %>%
dplyr::select('ref.date', 'ticker','ret.closing.prices') %>%
dplyr::filter(ticker == index_string) %>%
drop_na('ret.closing.prices')
rang = range(
index_ret$ret.closing.prices,
ticker_ret$ret.closing.prices
)
tib_returns <- ticker_ret %>%
inner_join(index_ret, by = 'ref.date')
print(tib_returns)
plot(
tib_returns$ret.closing.prices.y,
tib_returns$ret.closing.prices.x,
xlab = paste0("Market return", index_string),
ylab=ticker_string,
xlim=rang,
ylim=rang
)
capm_ticker <- lm(tib_returns$ret.closing.prices.x ~ tib_returns$ret.closing.prices.y)
abline(
capm_ticker$coef,
col = 2,
lwd = 3
)
title(
paste(
ticker_string,
"=",
round(capm_ticker$coef[1],4)," + ",
round(capm_ticker$coef[2],4),
"*SP500",
sep=""
)
)
fit <- dlmMLE(
y = tib_returns$ret.closing.prices.x,
parm = c(1, 1, 1),
x.mat = tib_returns$ret.closing.prices.y,
build = my_dlm,
hessian = T
)
se = sqrt(exp(fit$par))
mod_std = my_dlm(fit$par, tib_returns$ret.closing.prices.y)
mod_filt = dlmFilter(tib_returns$ret.closing.prices.x, mod_std)
mod_smot = dlmSmooth(mod_filt)
# colocamos o codigo original (com date -1) e não funcionou
# colocamos sem e parece estar tudo ok
date = tib_returns$ref.date
print(paste("isso", length(date)))
print(paste("esse", length(mod_filt$m[,1][-1])))
plot(
date,
mod_filt$m[,1][-1],
xlab = "day",
ylab = expression(alpha[t]),
type = "l",
main = ""
)
lines(
date,
mod_smot$s[,1][-1],
col = 2
)
abline(
h = capm_ticker$coef[1],
col=3
)
abline(h = 1, lty = 2)
plot(date,
mod_filt$m[,2][-1],
xlab = "day",
ylab = expression(beta[t]),
type = "l",
main = ""
)
lines(date,
mod_smot$s[,2][-1],
col=2
)
abline(
h=capm_ticker$coef[2],
col=3
)
abline(h = 1, lty = 2)
}
dlm_routine(asset_returns, ticker_string = "ITUB3.SA", index_string = "^BVSP")
evaluate_garch <- function(ticker, data, mean.model, variance.model) {
print(paste("Ticker", ticker, "usando os modelos"))
print(mean.model)
print(variance.model)
garch_spec <- ugarchspec(
mean.model = mean.model,
variance.model = variance.model,
distribution.model = "std"
)
garch_fit <- ugarchfit(spec = garch_spec, data = data$returns[-1])
print(paste("Visualizando erros do modelo para o ticker", ticker))
print(garch_fit@fit$coef)
tibble(ticker = ticker, fit = c(garch_fit), AIC=list(as.vector(infocriteria(garch_fit))))
}
garch_forecast <- function(garch_fit) {
garch_forecast <- ugarchforecast(garch_fit, n.ahead = 1)
garch_forecast %>%
plot
}
make_asset_garchs <- function(asset_returns, garchOrders) {
if (nrow(garchOrders) > 0 & nrow(asset_returns) > 0) {
tickers <- asset_returns$ticker %>%
unique
models <- tibble(ticker = NA, fit = NA, AIC = NA)
for (t in tickers) {
print(paste("Rodando garch para a ação ", t))
ret <- (asset_returns %>%
filter(ticker == t) %>%
drop_na('ret.closing.prices') %>%
mutate(
returns = ret.closing.prices
))
mean.model <- (garchOrders %>%
filter(ticker == t))$mean.model
variance.model <- (garchOrders %>%
filter(ticker == t))$variance.model
fit_result <- evaluate_garch(t, ret, mean.model, variance.model)
models <- models %>%
add_row(fit_result)
}
models
} else {
print("Assets e ordens dos garchs são necessários")
}
}
garchOrders <- tibble(
ticker = asset_returns$ticker %>%
unique(),
mean.model = list(
c(0, 0),
c(0, 0),
c(0, 0),
c(0, 0),
c(0, 0),
c(0, 0),
c(0, 0)
),
variance.model = list(
c(1, 2),
c(1, 2),
c(1, 2),
c(1, 2),
c(1, 2),
c(1, 2),
c(1, 2)
)
)
result <- make_asset_garchs(asset_returns, garchOrders)
for (t in unique(asset_returns$ticker)) {
print(paste("Info criteria do ticker", t))
print(infocriteria((result %>%
filter(ticker == t))$fit[[1]]))
}